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n IMAGE PROCESSING METHOD AND MEDICAL VIEWING SYSTEM FOR CHARACTERIZING 

SURFACES IN MEDICAL IMAGING" 

Description 

Field of the Invention 

The present Invention relates to an image processing method of characterizing 
boundaries of regions in an image for matching surfaces In medical imaging. The Invention 
further relates to a medical viewing system having means for carrying out this method and to a 
medical examination apparatus having such a medical viewing system. The Invention also 
relates to a computer program product having instructions for carrying out the method steps. 
The invention finds its application in the field of medical imaging and, more especially, In the 
field of x-ray medical imaging. 

When analyzing 3D medical images, It is often useful to match a discrete surface model 
to the boundary of an anatomical object. This can help, for example, in extracting the object 
from other unrelated objects for improving the visualization of said object Among many other 
applications, such discrete surface model allows to make clinical measurements such as linear 
dimensions, boundary surface area, cross-sectional area, enclosed volume of an organ. It Is 
also often desirable to compare or match active surface models with respect to one another. 

Background of the Invention 
A publication entitled "Edge Dipole and Edge Field for Boundary detection", pp.179-189, 
in "Part of the SPIE Conference on Hybrid Image and Signal Processing VI, ORLANDO, Florida, 
April 1998, SPIE Vol. 3389", by Toshiro KUBOTA eta/., proposes a framework for treating edges 
as directional dlpoles that induce a field around themselves. This paper treats edges not only as 
vectors but also as dlpoles, which Induce a vector field around them. Thus, said dipoles are 
called edge dipoles and the field is called edge field. An analogy is made between this concept 
and the Interaction of magnetic dlpoles with the magnetic field. The field induced by an edge 
dipole exhibits smooth continuous flows from toe positive pole to the negative pole. The field In 
turn can be used to rotate other dlpoles to align with the field. The approach of this paper 
provides an interaction mechanism between the dipole and the field, which has an influence on 
the edges. Also, the edge dipole field has circular flows diverging at the positive pole and 
converging at he negative pole, which provides a convenient mechanism for lateral Inhibition 
and discourages thick boundaries from forming. The dipoles interact with the field and align 
themselves into a smooth contour configuration. This concept is used in edge linking In order to 
complete missing image information. So, according to the disclosed method, it is needed to 
compute the influence of the dipoles within the region surrounding them in the images 
themselves. 

Summary of the Invention 

The present Invention has for an object to provide an image processing method of 
characterizing the boundaries of regions using edge dipoles without computation related to the 
image region Itself. The present Invention has further for an object to provide such an Image 
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processing method for analyzing an anatomical surface of Interest for comparing, matching or 
adapting a geometric model with anatomical object boundaries observable in a medical image, 
and/or a geometric model with another geometric model, and/or real image data with other 
image data of anatomical object boundaries observable in a medical image. 

Using the proposed method, substantial portions of the objects can be matched 
together in a patch wise approach. The present method is valid whether the objects are virtual, 
such as active models, or are real such as represented by real image date. A basic principle of 
the invention Is to represent a boundary of an object by dipoles alike a set of small elemental 
magnets. Each dipole characterizes an element of the boundaiy. Each dfpole has a strength 
proportional to the size of the boundary element, Is normal to said boundary element and 
oriented toward a direction regarded as outward direction with respect to the object limited by 

the boundary. 

The main steps of the method are recited In Claim 1. 

A medical viewing system for carrying out this method is further proposed. The medical 
viewing system can be Implemented as a specially programmed general-purpose computer. The 
medical viewing system can be a workstation. The present Invention yet further provides a 
computer program product having a set of instructions, when in use on a general-purpose 
computer, to cause the computer to perform the steps of the above-described method. The 
present invention still further provides a medical examination apparatus incorporating medical 
imaging apparatus, data processing system putting into practice the above-described method to 
process medical image data obtained by the imaging apparatus, and means for visualizing the 
image data produced by the method. The visualization means typically consists of a monitor 
connected to the data processing apparatus. Advantageously, the workstation and medical 
imaging system of the present Invention are Interactive, allowing the user to Influence clinical 
data that are evaluated and/or the manner In which evaluated data is to be visualized. 

Brief Description of the Drawings 
The invention and additional features, which may be optionally used to implement the 
invention to advantage, are described hereafter with reference to the schematic figures, where: 

Fig.lA illustrates dipoles associated to a boundary of an object and FIG.1B Illustrates 
the definition of the polynomial function parameters; FIG.1C illustrates a Dirac function; 

RG.2 Is a block diagram illustrating elements of a medical examination apparatus, 
incorporating a medical viewing system; 

Rg.3A Is a flow diagram showing the main steps of an image data processing method 
for characterizing a boundary and FIG.3B Is a flow diagram showing the main steps of an Image 
data processing method for matching two boundaries; 

FIG.4A illustrates a starting mesh model and RG.4B illustrates a step of fitting the 

model to a heart cavity; 

F1G.5 illustrates the mesh model that matches the heart cavity. 
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Description of Embodiments 

The invention relates to an image processing method for characterizing boundaries of 
regions in an image. The invention further relates to the application of this method to process 
medical images for improving the visualization of an anatomical surface of Interest The present 
Image processing method may be used for analyzing an anatomical surface of Interest for 
comparing, matching or adapting a geometric model with anatomical object boundaries 
observable in a medical Image, and/or a geometric model with another geometric model, 
and/or real image data with other image data of anatomical object boundaries observable In a 
medical image. Matching together substantial portions of the objects can be performed In a 
patch wise approach. The present Invention preferably makes, use of mathematical tools based 
on a polynomial function such as the Hermite transform, also known as Hermite wavelets. This 
mathematical tool Is succinctly described hereafter to facilitate computation steps that are 
advantageous to carry out the method of the invention. 

In a 2D or 3D Image, the Image Intensity f is defined in funcuon of the location of 
each image pixel or image voxel. So, in a 2D image, the scalar value image intensity f is given 
as function of the real positional coordinates x,y . A positive valued function w(x,y) is 
Introduced to define a fuzzy observation window W in the 2-D image. Geometric moments are 
defined as a set of quantities: M^,, • J"Jw(x, y) f(x, y) x m y" dx dy (l) 

The whole family of moments provides an alternative representation of the image flfoy) that 
fully and uniquely characterizes the intensity fi(x,y) within the observation window w(x,y) . In 
practice, lower orders of the data (m,n) are the most robust whereas higher orders are more 
delicate to compute numerically and are more noise-sensitive. These moments are powerful in 
coding morphological Information, but they have two shortcomings: 

a) When comparing two images together, within a similar observation window, it is not 
obvious how similarity/ dissimilarity measures of their moment representation relate to 
similarity/ dissimilarity of the original images. 

b) It is tedious to work out a comprehensive set of related measures that are Invariant 
under a transformation group. 

In order to solve these shortcomings a) and b), instead of using monomials x^y" in (1) 

one can use polynomials P^ a (jf, y) , which are mutually orthogonal within the observation 

window w(x,y) in the sense that: 

C nvA,.A, rf =* JJw(x,y)P nyi (x,^)P Irf ^(x,y)dx.dy (2) 
where each and are equal to 1 if the two Integers indices (m,m*) or (n,n') 
coincide and to 0 otherwise, and where > 0 is referred to as the norm of the polynomial 
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P „(*, y) . In feet It Is possible to build a complete set of orthogonal polynomials provided the 
wLow function w(x,y) Is never negative. Such a polynomial basis can provide least mean 
square approximation of an arbitrary function fl(x,y) in the form: 

sx^EWPn^yyc^ (3) 

ma 

where the 2-D Hermite coefficients are given by: 

f„ n = Jfw(x,y).P nwl (x,y).f(x,y).dxdy C 4a > 
Here, the coefficients f„, n replace the geometric moments defined in the previous 
paragraph. Those new coefficients are referred to as orthogonal moments. An Important 
property answers the requirement a) stated at the beginning of the current section. This 
theorem, called Parseval's theorem, says that the scalar products are conserved. So, 
considering any two real-valued Images f(x,y) and g(x,y) with respective orthogonal 
moments and gnMl , it is seen that taking the scalar products in transformed space, where 
the orders m and n play the role of coordinates, is the same as taking the product in x,y space: 
S^/C^ = Jfw(*,^x,y)g(x,y)dxdy (5) 

This relationship can be deduced from the orthogonality of the polynomial functions P^ . It is 
noteworthy that the normalization factors 1/C^ play In the m,n space, the role that the 

window function w(x,y) plays in the x,y space. 

The second shortcoming b) of the geometric moments Is that the moments do not 
transform easily with changes of referential. Therefore, it is required to simplify transformations 
of the orthogonal moments with change of referential. 

An advantageous way of doing this is to choose Hermite polynomials. This corresponds 
to taking for the window an isotropic Gaussian function centered about any reference 
point (£,*7)/ so that: 

W(X ' y)= ^- eXP( " 2p } 
This resulting transform Is called the Hermite transform. The corresponding 2-D polynomial 
basis can be expressed In term of products: 

P^y^H^JEUy) ™ 
where H m and H n are Hermite polynomials of orders m and n respectively. The corresponding 

(8) 

normalization factors are: = 2 m+n mln! 

Since an object of the invention is to compare images together, It is needed to easily 
deal with geometric transformations. 
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Dealing with translations Ss first described: The Hermfte coefficients f m ^ of function fQ 
can be expressed In terms of partial derivatives of a smoothed function F() obtained by 
convolving f() with the Gaussian kernel w(). Those skilled in the art can easily demonstrate 
that: 



/ CT ,=2 * .(-cr) 



dm+o e» 



dx a dy K 



This form allows to obtain the effect of translation by a simple application of Taylor 
theorem where the coefficients are directly related to the Hermite transform coefficients f^ n . 
The fact that a Gaussian smoothing was obtained to evaluate F implies that it is enough to use 
very few terms for toe Taylor expansions provided one stays within a disk of radius 1 or 2tr 
from the original window centre. 

Dealing with rotations is then described: To deaf with 2D rotations, the most suitable 
related transform Is that of Gauss-Laguerre transform described in the following publication by 
G. Jacovitrj and A. Neri; "Multiresolution Circular Harmonic Decomposition", IEEE Transactions 
on Signal Processing; vol. 48, pp. 3242-3247 (2000). This transform may be related to that of 
Hermite since the corresponding polynomials are related to solutions of the Quantum harmonic 
oscillator as described in the following publication by J J. Koenderink and A. van Doom; 
"Generic Neighbourhood Operators", IEEE Transaction on Pattern Analysis and Machine 
Intelligence, vol. PAMI-14, pp. 597-605 (1992). Hermite polynomials result from seeking 
solutions using Cartesian coordinates whereas the Laguerre polynomials appearing in the 
Gauss-Laguerre transform result from the use of polar coordinates. The essential thing one 
needs to know Is that, for any given order (pam+n) the two sets of solutions can be 
transformed Into each other. The Cartesian set of polynomials can be written as products 
Pm*(x,y)= a H B1 (x).H„(y) whereas the polar form (Gauss-Laguerre) are In the form 
fcos(*0)l 

e '* (r) tsto(*0)J < 10 > 

where Qp^r) is a function of the radial coordinate r = <Jx z + y 2 and 9 Is the angular polar 
coordinate. The order k of the circular harmonic factor spans all integers in the interval [0, p] 
with same parity as p (for k=0 only one circular harmonic exists, namely the constant 1). It Is 
easily seen that, for both representations, there are exactly (p+1) independent functions. In 
order to find the (p+1) 2 matrix that transforms the Hermite Into Gauss-Laguerre functions, it is 
enough to compare angular behaviours for very large r values. Elementary trigonometry can be 
used to express Pm. n (x,y) <* cos m (6).sin n (e) In terms of the circular harmonic functions referred 
to earlier. It Is also easily established that the same matrix can be used to transform 
coefficients from one basis to the other. One can also deduce the new normalization factors 
from this matrix because both set of functions form orthogonal bases. In order to revert from 
Gauss-Laguerre back to Hermite, the Inverse matrix Is used. 
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Rotation with an angle AO Is easily implemented because the coefficients Cp, k and S Pfk 
related to the Gauss-Laguerre functions Q p<k (r).cos(ke) and Q P(k (r).sin(ke) behave exactly like x 
and y components of a 2-D vector rotating with an angle kAO (for non-zero angular harmonic 
order k). Thus, a rotation would transform C^and Sp* to C p , k and S'pjsuch that 

( C "A = f cos < M<? > sin(M0) Y C pJk \ 

Which is totally equivalent to the complex variables equation: 

(C p * +^)=exp(^0)(C^ + «?„) < 12 > 
where I stands for the imaginary unit and AO stands for the rotation angle of the referential. 

Dealing with scale changes is further described: Changing scale implies changing the 
Gaussian window function. However, one can alternatively imply the scale changes in the 
approximating polynomials P^Q. Although this approach is relatively simple, It Implies using a 
least-mean square approximation Obtained within a given window for extrapolation Into a larger 
window or for interpolation within a smaller window. Scale changes are most easily operated on 
monomials. The approach proposed is just to express the polynomials as a sum of monomials, 
operating scale change and then reverting back to the orthogonal representation. For Hermife 











<\ 


0 


0 


0> 


H t (x) 
H 2 (x) 




H 2 (Xx) 


A = 


0 
0 


X 
0 


0 

X 2 


0 
0 






K ... j 






0 


0 


"V 



(13) 



and where A Is the lower triangular matrix allowing to express the column vector of Hermlte 
polynomials in terms of the column vector, of the powers of x. This allows to express H„(*x) as a 

20 linear combination of original Hermite polynomials H„(x) of equal or smaller order (and of same 
parity). The matrix A A A 1 can be computed once for all, it Is a sparse lower triangular matrix in 
which all non-zero coefficients are polynomials In X. of orders never exceeding the order of the 
line. Since 2-D Hermite polynomials are just products of ID polynomials, this scaling operation 
is dimensionally separable. Exactly the same transformation applies to Hermite coefficients. 

25 Because of the sparseness of the matrices involved and of dimensional separability, scale 
changes are not computationally costly. 

The above-described operations can be applied in images of more than two- 
dimensions: Most of the above-described operations easily generalize to any number of Image 
dimensions. The only basic change is how one deals with rotations. Spherical harmonics have to 

30 replace circular harmonics. For 3-D, the corresponding spherical harmonic formulae can be 
found in classical mathematics handbooks and the conversion matrices from Cartesian 
(Hermite) to polar form (Laguerre) can recursively be evaluated as is done In 2-D. However, a 
rotation of the referential implies the use of difficult spherical harmonics addition formulae. 
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Another formulation has been disclosed by the Clifford Algebra group of Ghent University by R. 
Delanghe, F. Sommen, V. Soucek, F. Brackx in "Clifford Algebra and Splnor-Valued Funcuons, A 
Function Theory of the Dirac Operator", Kluwer Academic Publishers (1992). For the 3-D case! 
the approach proposed by this publication just needs the introduction of unit called 
"quaternions" to Implement 3-D rotations of the Gauss-Laguerre coefficients, thus providing a 
generalization of the complex rotation factor exp(ikA0) as discussed In 2D case. Now, 
quaternions are entering the engineering field and are not considered to be exotic anymore in 
the Image processing community, as described by B, Horn in "dosed Form Solutions of 
Absolute Orientation using Unit Quaternions" in Journal of the Optical Society, vol. A4, pp. 639- 
642 (April, 1987). The basic advantage in using the above set of proposals Is to transform a set 
of image processing operations Into hopefully simple algebraic forms. These mathematical tools 
may be used by those skilled in the art In order to perform the computation steps of the image 
processing method of the invention. 

A method that applies this technique of characterizing a boundary using the function of 
Intensity with the Hermite coefficients is described hereafter. The present image processing 
method has steps for computing Hermite coefficients in order to characterize a curve line or a 
boundary of an object In an image window by the Hermite transform. This method permits of 
modelling and extracting Interfaces or boundaries of objects In medical Images. For performing 
this operation, it is possible to compute Hermite Transform Coefficients of an object boundary 
regarded as an oriented interface by sampling the corresponding boundary. There is no need to 
get region data as described in the cited "KUBOTA" publication. 

FIG.1A represents schematically a region of interest ROI of an object In a two- 
dimensional (2-D) image. The region of Interest ROI shows a 2-D boundary. A window W 
defines a portion SOI of the boundary that is considered. The 2-D boundary SOI is decomposed 
into small adjacent elementary linear elements, called segments, such as SI, S2, S3, S4, S5, 
respectively related to a predetermined number of boundary pixels, called reference pixels. ' 
Every segment such as S1-S5 of the boundary SOI in a 2D image Is represented by a respective 
dipole Dl, D2, D3, D4, D5. The arrows Dl to D5 illustrate the boundary dipoles. In 3-D images, 
the region of Interest ROI shows a 3-D surface of Interest SOI, called 3-D boundary, which is 
decomposed into small portions in a window W. Every small portion Is represented by a dipole 
of strength proportional to the area portion of the 3-D boundary and oriented along the 
outward normal direction, at a defined centre of the small portion. 

From such a representation, the Hermite transform can be computed within the window 
W as schematically represented in FIG.1A for a 2-D Image. The intensity of the boundary SOI in 
the window W is defined by a function of intensity f(x, y) where x,y are the coordinates of 
the reference pixels at the center of the segments or small portions. Assuming that the window 
centre is at the origin and that distances are normalized so that <rl2=l, the Hermite transform 
coefficients related to a dipole vector (lx,Iy) at point (x,y) can be expressed as: 
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f m n - w(x,y).(Ix.Hm+l(x).Hn(y) + ¥"m (x).H n+ l(y)) W 
The Hermite transform of the interface dipoie layer can be approximated as the sum of such 
contributions over all discrete interface elements in the window. Generalization to more than 2- 
D is again straightforward. L _ c 

A first data, denoted by n, called spatial resolution, is related to the number of said 
segments in which the boundary is decomposed in said window W. A second date, denoted by 
m, called intensity resolution/ is related to the number of possibie gray ^ds tor tee pbcds o 
the boundary In said window W. The definition of the data n, m is illustrated by R6.1B. Near 
point O, the spatial resolution n is low, meaning the number of segments Is small while the 
segments are long; and the intensity resolution is low, meaning the number of gray levels «s 
small, instead near point p, tee spatial resolution n is high, meaning the number of segments * 
large while the segments are small; and the intensity resolution is high, meaning the number of 
gray levels is large. The date n, m may be chosen a priori* the user. So, using the values of 
n m near point O, the boundary SOI In the window W Is characterized by a few number of long 
segments and a few number of gray levels. Instead, using the values of n, m near point p, the 
boundary SOI in the window W is characterized by a large number of small segments and a 

large number of gray levels. 

In the present method, instead of regarding a boundary pixel directly, said boundary 
pixel is first coded using the above described technique. Each segment Is characterized by a 
dipoie, such as Dl, D2, D3, D4, D5, located at the centre x,y of the segment, oriented along 
the outward normal direction of the object and having a strength that Is defined using the 
spatial and intensity data related to tee segment. The coding data Include the first information 
related to the spatial resolution n and the second Information related to the Intensity resolution 
,ri the window W of the region of interest ROI. The spatial resolution n is considered as a first 
dimension and tee intensity resolution m Is considered as a second dimension for determining 
coding data related to a pixel at the location x,y . Alternately as shown in HG.1C, the segment 
may be coded by a Dirac Intensity function having a highly positive value, outside the object, 
with respect to an axis directed along the dipoie, and a highly negative value Inside tee object. 

The major difference between the present invention and tee technique disclosed in tee 
publication "KUBOTA etat cited above, Is teat instead of broadcasting the effect of each d.pole 
,n a portion of tee image as was done in the known technique, the effect of each such dipoie »s 
integrated into transformed coefficients. The effects of those dipoles, such as D1-D5, w,th,n tee 
window W, are cumulated as coefficients of an image transform procedure. The polynomial 
Hermite transform Is advantageously used. As above-described, said Hermite transform .s 
, capable of coding the boundary Inside tee window using a finite number of coefficients. The 
contribution of each dipoie element can be computed analytically knowing Its position within the 
window, Its strength and its orientation. The coding date form a series of two-dimensional (2-D) 
coefficients for defining said polynomial function. The number of 2-D coefficients is small when 
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the coding data are chosen near point O of FIG.1B, while the number of 2-D coefficients is large 
when the coding data are chosen near point p of RG.1B. The number of 2-D coefficients may 
be chosen within a range in which It may vary slowly, since the coding data n, m may be each 
chosen In a range In which they may vary slowly. Accoiding to the method of the invention, the 
5 boundary Is characterized by said polynomial function with said 2-D coefficients in the window 
W. If the polynomial function is constructed with a great number of such 2-D coefficients, the 
representation of the boundary is smooth. If the polynomial function is constructed with a small 
number of coefficients, the representation of the boundary is coarse. Hence said polynomial 
function with said 2-D coefficients permits of classifying the representation of the boundary in a 
10 hierarchical manner, from coarse to fine, according to the number of coefficients used to 
construct the polynomial function. The medical viewing system and an Image processing 
method of the present invention permit to minimize the load of computation produced by the 
KUBOTA eta/, approach. 

FIG.3A Is a flow diagram of main steps of the present method, comprising: 
In step SO, acquisition of image data: The image data input to the method can be for 
example 3-D computed tomography image data obtained for a subject heart. The medical 
image date consists of a large number of date relating to points, each corresponding to a 
respective position within the patlenfs body. The image data can be submitted to preprocessing 
in order to eliminate noise. The method further comprises steps of: 

In step SI, computation of image date of an object surface denoted by Bl or B2. In 
step SI, the outer surface of the heart musde, denoted by Bl, is identified from within the 
image date via a segmentation process as illustrated by the segmented 2-D curved line or 3-D 
surface SOI in Fig.lA. In another technique, the 3-D surface, denoted by B2, may be obtained 
as an active model providing a best fit to the heart muscle, or other anatomical object under 
consideration. Techniques for producing an active model of an anatomical object are well- 
known, for example by the description In the publication entitled "General Object 
Reconstruction Based on Simplex meshes" by Herve Delingette, In the International Journal of 
Computer Vision, 32, 111-142, 1999. Such a model Is Illustrated by FIG.4A and FIG.4C. 

In step S2, decomposition of the boundary into 2-D elements of a 2-D curve line or 
small portions of a 3D boundary, such as Bl or B2, as above described In reference to FIG.1A. 

In step S3, characterization of the 2-D or 3-D boundary elements for example by 
dlpoles or Dirac functions as illustrated by FIGS.1A-1C; and computation of corresponding 
Hermlte coefficients according to formulae (4a) or (4b). 

In step S4, cumulating the effects of the dipoles Dl to D5 within the window W as 
coefficients of an image transform procedure. For example, one can use a Hermite transform 
that is capable of coding an Image Inside the window W using a finite number of coefficients. In 
this step, the contribution of each dipoie element Dl to D5 in the window W Is computed 
analytically knowing Its position within the window, Its strength and its orientation according to 
formulae (3), (4), (6) and (7). 
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In step S5, computation of the Hermite coefficients for translation, rotation and scale 
cange according to formulae (9), (11), (12), (13). 

The Hermite transform of a boundary can now be used for different purposes. The 
present invention proposes a method having the steps described above and having further 
steps for matching two surfaces. For example this method has further steps to compare and/or 
to match a discrete model interface to another discrete interface; or to compare and/or to 
match a discrete interface to a real image. It is noted that the above proposals can substantially 
Improve the performance and computational complexity of matching procedures. None of the 
above procedures requires to know a mapping from initial point position (prior to 
transformation) to a target position (after transformations) as is needed by the cited algorithms 
disclosed in the publication by B. Horn. So, the invention further relates to applying the first 
steps of characterizing the function of Intensity of a boundary with the Hermite coefficients for 
comparing or matching two objects defined by curve lines or surfaces in corresponding 
windows. The operation of comparing or matching the boundaries of said two objects can be 
done only using the respective Hermite coefficients relating to their boundaries. The 
transformation needed to match a portion of interface with another one (registration or motion 
estimation) can be estimated for a pair of objects of same nature such as model with model or 
real data with real data; or of dissimilar nature such as model to data or data to model. The 
proposed method can be applied to medical image processing in order to enabling improved 
matching of surfaces such as an active model with anatomical object boundaries observable in 
a medical image, and/or an active model with another active model, and/or real image data 
with other image data of anatomical object boundaries observable in a medical image. 

Hence, this method comprises steps for computing Hermite coefficients in order to 
diaracterize.corresponding boundaries of two objects- in an image window by the Hermite 
transform, including Hermite coefficients for translation, rotation, and scale changing and 
further steps for comparing or matching said two boundaries. Here, the template (or block-) 
matching, which Is well known to those skilled In the art, can be transformed Into matching of 
transformed coefficients. Rather than making exhaustive searches for translations, rotations, 
and scale changes in said template (or block-) matching operations, the mean square matching 
error can be written as an algebraic function of the transformation parameters. Optimal match 
between two boundaries can then be found by conventional gradient descent Furthermore, it 
can be hoped that the simple algebraic form of the matching error can be exploited to 
efficiently locate all the local minima. A rather interesting feature of the Hermite transform is 
that the set of coefficients for which one of the indices is 0 corresponds to the Hermite 
Transform of a "Weighted Radon" projection along axes corresponding to the 0 index. This may 
be useful, for example, for matching a 3-D pattern with a 2-D projection pattern. 
Referring to Fig.3B that is a flow diagram, this method comprises: 
In steps Tl, T2: respectively acquiring image data of first and second Images. In the 
example described hereafter, in the first image the object is an organ represented by real data, 
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and in the second image, the object Is a virtual model of the organ having a boundary to be 
compared with a boundary of said organ. The interface of the virtual model is a binary region. 

In steps T3, T4: respectively segmenting said images to provide Image data of the 
corresponding boundaries respectively denoted by Bl, B2 in the first and second images. 
5 In step T5: calculating Hermite coefficients of the Hermite transforms for boundaries Bl 

and B2, as described in reference with steps S2, S3, S4 illustrated by FI6.3A. 

In step T6: deducing, as in step S5 illustrated by RG.3A, coefficients resulting from 
rotation, scale-change or translation of one object with respect to the other. The computations 
only require knowledge of the geometric transformation and of the Hermite coefficients before 
10 transformation. 

In step T7: estimating the transformation needed to match the portion of boundary or 
interface of the model with the portion of boundary or Interface of the organ, only using their 
respective Hermite coefficients. 

Another application Is for performing registration or for motion estimation. Another 
15 application consists in smoothing an interface model through adapting It to a flat, or another 
appropriately shape target interface. RGs.4A illustrates an active model to be matched to heart 
cavity as shown In FIG.4B. This model is deformed using Internal and external forces In such a 
way as to fit the heart cavity. In order to determine the best match, the model surface 
boundary, as shown In R6.5, is compared to a segmented surface boundary of the heart cavity 
20 using the above-described method. In this first example, the above-described steps are used 
for estimating the transformation needed to match a portion of interface of objects of dissimilar 
nature, such as a model compared to real data or such as real data compared to a model. The 
above-described steps can also be used for estimating the transformation needed to match a 
portion of interface of a pair of objects of same nature, such as a model compared to a model 
25 or such as real data compared with real data. 

Then, each reference point of the reference surface Is processed in turn and the normal 
to the reference surface Is calculated at each reference point. All operations Involved can be 
done In computationally efficient manner. The basic technique proposed here is very generic; 
more than the cited operations may be carried out. 

In general, the described method above, which is based on the technique of 
characterizing the function of intensity of a boundary with the Hermite coefficients of the 
Hermite transform of the boundary, can be used for purposes such as: 

Adapting one contour to an Image, which would consist in seeking a maximum of the 
correlation. This operation adapts the position of the dipoles in such a way as to maximize 
correlation of an interface with another Interface called target interface. This can be used for 
example to model large portions of the boundary of an anatomical object using an active 
contour or active surface as model. This can be done by changing the dipole positions within 
the window (free-form adaptation) or by seeking optimal rotation, scale change and/or 
translation (global adaptation within the window). In fact, the operation of comparing or 
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correlating the boundaries of said two objects can be done only using the respective Hermite 
coefficients relating to their boundaries. 

Smoothing an interface model: this is possible through adapting it to a flat or 
appropriately shaped target interface. Smoothing such a boundary is performed by minimising 
the amplitude of higher order terms of the Hermite transform. 

In another application, this method can be used to study noise or texture in images. 
This operation can be carried out using estimation of their autocorrelation. This consists in 
studying how the Image correlates to a translated version of itself. An extension of this idea is 
to characterize the way the image correlates to itself when it undergoes, translation, scale 
changes and rotations. One should therefore estimate: 

R (py (f(T~ 1 x)f(Tx}j , where T and T 1 are a transformation and its inverse. The 

windowed correlation product can be expressed algebraically in terms of the parameters 
defining the group of the geometric transformations T. If this group includes rotation, scaling 
and translation, the resulting autocorrelation function Is rotation, scale and translation 
independent If this is expressed in a functional form, coefficients characterizing this functional 
form can be extracted. For rectangular windows and when limiting T to mere translations, the 
functional form of the autocorrelation is characterized by the power spectra, distribution (square 
of moduli of the Fourier transform) through the Wiener Khintchine theorem as described by 
W H Press, SA. Teukolsky, W.T. Vetterling and B.P. Flannery; "Numerical Recipes in C, The Art 
of Scientific Computing", Cambridge University Press, 1992 (2« d Edition). Such characterisations 
of autocorrelation functions can be used to characterizing texture or the spectral behaviour of 
noise, and or providing an Invariant set of characteristics of shapes for pattern recognition. 

The above method for coding boundariescan be extended to other structures and to 
many applications. For example, oriented linear or tubular structures can be coded as a set of 
quadripoles. It is possible to put this transform in one of two forms: Hermite Cartesian form 
that is best suited to deal with translation or scale-change, and/or Gauss-Laguerre Polar form 
which is best suited to deal with rotations. Conversion from one form to the other is not 
computationally costly if one concentrates on low orders (/.e. when neglecting high frequency 
noise-prone terms in the expansions). Implementing these basic transformation operations can 
be done in a computationally efficient manner and can also be put in the form of compact 
algebraic expressions. Many practical applications can be envisaged as outlined in the previous 
section for putting the 3D rotation transformations Into a mathematically correct framework and 
or characterizing the Generalized Autocorrelation function and relating it with invariant theories 
used in pattern recognition. An important advantage of Hermite transform is the ease and 
efficiency with which it can be put In forms that are particularly suited to deal with translation, 
rotation and scale change. Once the coefficients of the Hermite transform are known, one can 
deduce the coefficients resulting from rotation, scale-change or translation of the object The 
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computations only require knowledge of the geometric transformation and of the Hermite 
coefficients before transformation. 

Ffg.2 shows the basic components of an embodiment of an Image viewing system in 
accordance to the present invention, incorporated In a medical examination apparatus. As 
Indicated schematically In Rg.2, the medical examination apparatus typically includes a bed 10 
on which the patient lies or another element for localizing the patient relative to the imaging 
apparatus. The medical imaging apparatus may be a CT scanner 20. The image data produced 
by the CT scanner 20 is fed to data processing means 30, such as a general-purpose computer, 
that carries out the steps of the method. The data processing means 30 is typically associated 
with a visualization device, such as a monitor 40, and an input device 50, such as a keyboard, 
pointing device, etc. operative by the user so that he can interact with the system. The 
elements 10-50 constitute a medical examination apparatus according to the invention. The 
elements 30-50 constitute a medical viewing system according to the invention. The data 
processing device 30 is programmed to Implement a method of processing medical image data 
according to invention. In particular, the data processing device 30 has computing means and 
memory means to perform the steps of the method. A computer program product having pre- 
programmed instructions to carry out the method may also be Implemented. 

The present invention Is applicable regardless of the medical imaging technology that is 
used to generate the initial data. For example, when seeking to visualize the heart magnetic 
resonance (MR) coronary angiography may be used to generate 3D medical image data in a 
non-invasive manner. See, for example, "Non-invasive Coronary Angiography by Contrast- 
Enhanced Electron Beam Computed Tomography" by Achenbach etal, in Clinical Cardiology, 21, 
323-330, 1998. The Achenbach et al article includes useful information regarding optional data 
processing steps that can be applied to the medical image data, for example, segmentation to 
enable a representation of certain anatomical features in isolation from others, etc These 
steps can be applied in the method of the present invention. The present invention Is applicable 
regardless of the way in which a surface of interest is modeled, whether via use of a reference 
simplex mesh, or In some other way. Various modifications can be made to the order in which 
processing steps are performed in the above-described specific embodiment. The above- 
described processing steps applied to medical image data can advantageously be combined 
with various other known processing/visualization techniques. The drawings and their 
description hereinbefore Illustrate rather than limit the invention. It will be evident that there 
are numerous alternatives that fall within the scope of the appended claims. Moreover, 
although the present Invention has been described In terms of generating image data for 
display, the present Invention is intended to cover substantially any form of visualization of the 
image data including, but not limited to, display on a display device, and printing. Any reference 
sign in a daim should not be construed as limiting the claim. 
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CLAIMS: 

1. A medfcal viewing system comprising data acquisition means for acquiring image data 
of an object in an image and processing means for characterizing a boundary of the object 
comprising computing means for identifying the object boundary and reference points of said 
object boundary within an observation window; and for decomposing said object boundary into 
boundary unit elements, centred at the reference points; 
this medical viewing system further comprising means for: 

coding each reference point using coding data including spatial information and 
Intensity information relating to said reference point and the corresponding unit element; 

computing transform coefficients from said coding data; 

representing said object boundary by a polynomial transform function using a finite 
number of said coefficients; and . 

visualizing object Images and/or processed images. 

2. The medical viewing system of claim 1, wherein the means for coding the reference 
points yield spatial the Information based on the size of the unit elements; and yield the 
Intensity information based on the number of intensity levels along the boundary. 

3. The medical viewing system of one of Claims 1 or 2, wherein the means for coding the 
reference points represent each boundary unit element by a dlpole of strength proportional to 
the size of the boundary unit element, centred at the reference points and oriented along the 
outward normal direction to the object boundary; and 

the means for yielding the transform coefficients compute said coefficients from said 
dipoles for further representing said object boundary by a polynomial transform function using 
a finite number of said coefficients. 

4. The medical viewing system of one of Claims 1 to 3, wherein the transform coefficients 
are Hermite transform coefficients. 

5. The medical viewing system of claim 4, wherein the computing means for computing 
transform coefficients from said dipoles comprise means for integrating the effect of each dlpole 
into the transform coefficients, the contribution of each dipole being computed analytically 
knowing Its position within the window, its strength and its orientation. 

6. The medical viewing system of one of Claims 1 to 5, comprising means for smoothing 
the boundary of an object by varying the number of used transform coefficients. 

7. The medical viewing system according to one of Claims 1 to 5, comprising means for 
comparing two corresponding boundaries of two different Images of an object, by comparing 
their transform coefficients. 

8. The medical viewing system of claim 7, comprising means for matching two 
corresponding boundaries of two different images of an object, by matching their transform 
coefficients. 
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9 The medical viewing system of claim 8, wherein the means for matehing the transform 
coefficients of said boundaries comprises means for applying translation, rotation, scale change 
to the transform coefficients of one of the corresponding boundary of the object in one .mage 

10 The medical viewing system of claim 8, comprising means for performing correlation of 
the corresponding boundaries of an object In different images by comparing the correspondmg 
transform coefficients of the boundaries in the images and matching said transform coefficients. 
U The medical viewing system of one of claims 6 to 10, comprising means for acquinng 
real image date for forming the two images of the object; or virtual Image date for forming the 
two images of the object; or real image date for forming the first image of the object and 
virtual image data forforming the second image of the object; and 

comprising means for comparing or matching or correlating corresponding boundanes 
of the different images of the object, by comparing or matehing or performing correlation of 
corresponding coefficients of the real or virtual boundaries, including when necessary means for 
applying translation, rotation, scale change to the transform coefficients of at least one of the 
15 corresponding boundary of the object in one image. 

11 An image processing method to cause the data processing means of the medical 
viewing system of one of Claims 1 to 11 to perform steps of acquiring image date of an object 
in an image and characterizing a boundary of the object, comprising identifying the object 
boundary and reference points of said object boundary within an observation window; and 
20 decomposing said object boundary into boundary unit elements, centred at the reference 
points; this method further comprising steps of: 

coding each reference point using coding date including spatial informabon and 
intensity informabon relating to said reference point and the corresponding unit element; 
computing transform coefficients from said coding date; 
25 representing said object boundary by a polynomial transform function using a finite 

number of said coefficients; and 

visualizing object images and/or processed images. 
12. A medical examination apparatus comprising acquisition means for acquiring medical 
image date, and a medical viewing system according to one of Claims 1 to 11. 
30 13 A computer program product having a set of instructions, when In use on a general- 
purpose computer, to cause the computer to perform the steps of tee method according to 
Claim 11. 
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" IMAGE PROCESSING METHOD AND MEDICAL VIEWING SYSTEM FOR CHARACTERIZING 
SURFACES DM MEDICAL IMAGING" 

Abstract; 

The Invention relates to a medical viewing system comprising data acquisition means 
for acquiring Image data of an object In an image and processing means for characterizing a 
boundary of the object. This system comprises computing means for Identifying the object 
boundary and reference points of said object boundary within an observation window; and for 
decomposing said object boundary into boundary unit elements, centred at the reference 
points. This medical viewing system further comprises means for coding each reference point 
using coding data Including spatial, information and intensity information relating to said 
reference point and the corresponding unit element; computing transform coefficients from said 
coding data. A finite number of said coefficients is used for representing said object boundary 
by a polynomial transform function. This system has means for comparing, matching, 
correlating, smoothing corresponding boundary portions of the object In different images, only 
by performing operations on the transform coefficients. These operations yield rotation, 
translation, scale change transform coefficients. Said images may be either virtual or real. 
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